#pragma rtGlobals=1		// Use modern global access method.
#include <WaveSelectorWidget> // Includes the WaveSelectorWidget Procedure

//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//
//																	 //
//																	 //
//  						GUI Procedures								 //
//																	 //
//																	 //
//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//



//--------------------------------------------------------------------------------------------------------------
// Main control panel recreation macro
//--------------------------------------------------------------------------------------------------------------

Window faMainPanel() : Panel
	PauseUpdate; Silent 1		// building window...
	NewPanel /N=faMainPanel/W=(455,55,1275,375) as "Factor Analysis"
	Button makeMatrix,pos={12,292},size={80,20},proc=ButtonProc,title="New Matrix"
	Button mainCancel,pos={756,292},size={50,20},proc=cancelButton,title="Cancel"
	SetVariable matrixSelected,pos={14,40},size={182,16},bodyWidth=150,title="Matrix"
	SetVariable matrixSelected,value= selectMatrix,live= 0, proc=setVarCntrl
	Button deCompButton,pos={12,74},size={90,20},proc=deCompButton,title="Decompose"
	Button compFact,pos={119,74},size={65,20},proc=factCompressButton,title="Compress"
	ListBox eigenValues,pos={205,42},size={600,235},listWave=root:eigenListBox
	ListBox eigenValues,selWave=root:eigenValueList,mode= 9, titleWave = root:ListBoxTitles, widths={90}
	CheckBox check0 title="Normalize",mode=0, variable=root:normalizeDat, pos={12,100}, proc=normRadioButton
	SetWindow faMainPanel hook(xButton) = exitButtonHook
EndMacro



//--------------------------------------------------------------------------------------------------------------
// Add Menu item
//--------------------------------------------------------------------------------------------------------------

Menu "Spectrometry"
	"Factor Analysis", factAnaInterface()
end



//--------------------------------------------------------------------------------------------------------------
// Start up procedure to build main control panel and do data folder and 
// variable initation
//--------------------------------------------------------------------------------------------------------------

Function factAnaInterface()
// Sets the data folder to root
setdatafolder root:

// Checks is the Factor_Analysis data folder exists. If it doesn't then create it
If(dataFolderExists("Factor_Analysis"))
	Variable/G root:normalizeDat = 0
else
	NewDataFolder Factor_Analysis
	Variable/G root:normalizeDat = 0
endif

// Waves needed for the main control panel. Will be used to display eigen values
// After decomposition
Make/O/T/N=0 defaultListBox
Make/O/T/N=(0,6) eigenListBox
Make/O/T/N=6 ListBoxTitles = {"Eigen Values", "RE", "IE", "XE", "IND","%SL"}
Make/O/B/N=0 eigenValueList=0

// Global variables needed for various procedures
String/G nameMatrix="Name your matrix", selectMatrix = "Matrix Name Goes Here", numbfact = "0;"

// Checks to see is the main panel is already open. If it is then bring it to the front
// if not then creates it
if(wintype("faMainPanel") != 0)
	DoWindow/F faMainPanel

elseif(wintype("faMainPanel") == 0)
	Execute "faMainPanel()"

endif

end



//--------------------------------------------------------------------------------------------------------------
// Procedure to add waves to the selected waves list
//--------------------------------------------------------------------------------------------------------------

Function updateSelected()	
// Variable declaration							
String selectedWaves = WS_SelectedObjectsList("dSelect", "WaveSelectorList"), waveListHolder = ""
Variable q=0

// do-while loop to remove the path information from the selected waves list for display
// in the selected waves list box
Do
	String selectedWavesHolder = StringFromList(q,selectedWaves, ";")[strsearch(StringFromList(q,selectedWaves,";"), ":" ,strsearch(StringFromList(q,selectedWaves, ";"), ":", inf, 1)-1, 1)+1,inf]
	waveListHolder = waveListHolder +  selectedWavesHolder + ";"
	q+=1
	
while(q<itemsinlist(selectedWaves))

// Duplicates the already selected waves wave for concatenation
Duplicate selectedWaves_Holder, firstWave

// Creates a wave to hold the current wave selections
Make/T/N=(ItemsInList(waveListHolder, ";")) waveHolder
waveHolder = StringFromList(p, waveListHolder, ";")

// Dupicates the alread selected wave paths for concatenation
Duplicate selectedWaves_wave, firstwave_path

// Creates a wave to hold the path to the waves in the current wave selection
Make/T/N=(ItemsInList(selectedWaves, ";")) selectWave
selectWave = StringfromList(p, selectedWaves, ";")

// Concatenates the waves to add the current selection to the previous selection
Concatenate/O/KILL/NP {firstWave, waveHolder}, selectedWaves_Holder
Concatenate/O/KILL/NP {firstWave_path, selectWave}, selectedWaves_wave

// Redimensions the wave that holds the selection status of items in the
// Selected waves list box to match the new selection size
Redimension/N=(numpnts(selectedWaves_Holder)) numSelectedwaves

// Clears the selection to make way for a new selection
WS_ClearSelection("dSelect", "WaveSelectorList")

end



//--------------------------------------------------------------------------------------------------------------
// Procedure for removing selected waves from the selected waves list box
//--------------------------------------------------------------------------------------------------------------

Function removeSelected()
// Wave and variable declaration
WAVE numSelectedWaves,selectedWaves_Holder,selectedWaves_wave
Variable holder=0, n=0

// do-while loop to delete the selected waves from the selected waves list box
// and the various waves
Do 
	holder = numSelectedWaves[n]

	If(holder >= 1)
		Deletepoints n, 1, numSelectedWaves
		Deletepoints n, 1, selectedWaves_Holder
		Deletepoints n, 1, selectedWaves_wave
		n=0
		continue
	else
	
	endif
	
	n+=1

While(n<numpnts(numSelectedWaves))

end



//--------------------------------------------------------------------------------------------------------------
// Procedure for selecting the waves to form the data matrix
//--------------------------------------------------------------------------------------------------------------

Function waveSelect()

// Sets the data folder to root
SetDataFolder root:

// Creates the data selection control panel
NewPanel/N=dSelect/W=(181,179,753,510) as "Data Selection"

// list box control doesn't have any attributes set on it for use by wave selector widget
ListBox WaveSelectorList,pos={9,13},size={273,241}

// This function does all the work of making the listbox control into a
// Wave Selector widget. Note the optional parameter that says what type of objects to
// display in the list. 
MakeListIntoWaveSelector("dSelect", "WaveSelectorList")

// String to hold the selected waves
string selectedWaves = WS_SelectedObjectsList("dSelect", "WaveSelectorList")

// Wave to hold the path to the selected waves
Make/T/N=(ItemsInList(selectedWaves, ";")) selectedWaves_wave

// Wave to hold the selected waves for display in the selected waves list box
Duplicate selectedWaves_wave, selectedWaves_Holder

// Wave to hold the selection status of the items in the selected waves list box
Make/O/B/N=(numpnts(selectedWaves_Holder)) numSelectedWaves=0

// Transfers the selected waves from the string to the wave
selectedWaves_wave = StringFromList(p, selectedWaves, ";")

// Creates the rest of the control panel needed for selecting the waves to create the Matrix
ListBox WaveSelectedList,pos={292,13},size={273,241},listwave=selectedWaves_Holder, mode=10, selWave=numSelectedWaves
Button cancel,pos={9,300},size={50,20},proc=cancelButton,title="Cancel"
Button select,pos={69,300},size={50,20},proc=selectButton,title="Select"
Button deselect,pos={129,300},size={60,20},proc=deselectButton,title="Remove"
Button matrix,pos={199,300},size={110,20},proc=makeMatrix,title="Create Matrix"
SetVariable matrixName,pos={319,300}, size={245,20}, title="Matrix Name", value=nameMatrix
SetWindow dSelect hook(xButton) = exitButtonHook

end


//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//
//																	 //
//																	 //
//  					   Control Panel Procedures						 //
//																	 //
//																	 //
//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//



//--------------------------------------------------------------------------------------------------------------
// Button control procedure to add to the selected waves
//--------------------------------------------------------------------------------------------------------------

Function selectButton(ctrlName) :ButtonControl
string ctrlname
updateSelected()
end

Function normRadioButton(ctrlName,checked) :CheckBoxControl
string ctrlname
variable checked
NVAR normalizeDat = root:normalizeDat
normalizeDat = checked
 
end

//--------------------------------------------------------------------------------------------------------------
// Button control to remove waves from the selected waves
//--------------------------------------------------------------------------------------------------------------

Function deselectButton(ctrlName) :ButtonControl
string ctrlname
removeSelected()
end



//--------------------------------------------------------------------------------------------------------------
// Button control to make the matrix from the selected waves
//--------------------------------------------------------------------------------------------------------------

Function makeMatrix(ctrlName) :ButtonControl
string ctrlname

// Variable declaration
SVAR nameMatrix
String wavepaths = "", fullNameMatrix = nameMatrix + "_Factor_Analysis"

// Sets the text wave to be used to update the displayed list
WAVE/T selectedWaves_wave

// Variable initilization for the do-while loop that updates the selected waves list
Variable n=0

// Error handling to detect if user pressed create matrix without selecting any waves
if(numpnts(selectedwaves_wave) == 0)
	PauseUpdate; Silent 1		// building window...
	NewPanel /N=selectionError/W=(150,77,406,167) as "Selection Error"
	TitleBox errorMessage,pos={17,28},size={344,16},title="Please select at least one data wave."
	TitleBox errorMessage,fSize=14,frame=0
	Button buttonOk,pos={16,57},size={50,20},title="Ok",proc=okButton
	SetWindow selectionError hook(xButton) = exitButtonHook
	PauseforUser selectionError
else
	// Error handling to check if the chosen matrix name is valid
	if(CheckName(nameMatrix, 6) != 0 || CheckName(fullnameMatrix, 11) != 0)
		PauseUpdate; Silent 1		// building window...
		NewPanel /N=nameError/W=(150,77,476,167) as "Illegal Matrix Name"
		TitleBox errorMessage,pos={17,28},size={344,16},title="Invalid Matrix name. Choose another Matrix name."
		TitleBox errorMessage,fSize=14,frame=0
		Button buttonOk,pos={16,57},size={50,20},title="Ok",proc=okButton
		SetWindow nameError hook(xButton) = exitButtonHook
		PauseforUser nameError
	else
		// Creates path to the new matrix
		string matrixPath = "root:Factor_Analysis:" + nameMatrix
		
		// Tests for existance of a matrix with the same name
		if(exists(matrixPath) == 0)
			// do-while loop that creates string list from the text wave to pass to the dataMatrix function
			Do
				wavepaths = wavepaths + selectedWaves_wave[n] + ";"
				n+=1
			while(n<numpnts(	selectedWaves_wave))
		      
		      // Calls the dataMatrix function with the inputed parameters and then kills 
		      //the select waves control panel
			dataMatrix(wavepaths, nameMatrix)
			cancelButton("cancel")
		
		// Error handling in case there is a name conflict with the user enter matrix name	
		else
			PauseUpdate; Silent 1		// building window...
			NewPanel /N=nameError/W=(150,77,476,167) as "Illegal Matrix Name"
			TitleBox errorMessage,pos={17,28},size={344,16},title="Matrix already exists. Choose another Matrix name."
			TitleBox errorMessage,fSize=14,frame=0
			Button buttonOk,pos={16,57},size={50,20},title="Ok",proc=okButton
			SetWindow nameError hook(xButton) = exitButtonHook
			PauseforUser nameError
		
		endif
	
	endif

endif

end



//--------------------------------------------------------------------------------------------------------------
// Button Control procedure to initate the wave selection control panel
//--------------------------------------------------------------------------------------------------------------

Function ButtonProc(ctrlName) : ButtonControl

String ctrlName
// Initates wave selection control panel and pauses to keep user from interacting
// with main control panel
waveselect()
PauseforUser dSelect

End



//--------------------------------------------------------------------------------------------------------------
// Button Control procedure to handle the cancle button on all control panels
//--------------------------------------------------------------------------------------------------------------

Function cancelButton(ctrlName) :ButtonControl
string ctrlName

// Checks for which window the cancel button was pressed in and then 
// kills that window and cleans up and waves or variables no longer needed
If(stringmatch("mainCancel", ctrlName))
	
	setDataFolder root:
	KillWindow faMainPanel
	//setDataFolder root:
	KillWaves/Z root:defaultListBox,root:eigenValueList,root:eigenListBox,root:ListBoxTitles
	KillVariables/Z V_flag
	KillStrings/Z nameMatrix, selectMatrix, numbfact, S_name
	
elseif(stringmatch("cancel", ctrlName))
	KillWaves/Z root:selectedWaves_wave, root:selectedWaves_Holder, root:numSelectedWaves
	if(wintype("dSelect") != 0)
		KillWindow dSelect
	endif

endif

end



//--------------------------------------------------------------------------------------------------------------
// Button Control procedure to handle the Decompose button
//--------------------------------------------------------------------------------------------------------------

Function deCompButton(ctrlName) :ButtonControl
String ctrlname

// Global variable to pass to the decomp procedure
SVAR selectMatrix = root:selectMatrix

// Sets the data folder to the Factor_Analysis folder were all raw matrices are stored
setdatafolder root:Factor_Analysis

// Error handling in case the entered matrix name doesn't exists. 
// In the future this should be replaced by another wave selection control panel
// that displays only the waves in the Factor_Analysis  data folder
if(exists(selectMatrix) == 0)
	PauseUpdate; Silent 1		// building window...
	NewPanel /N=nameError/W=(150,77,536,167) as "Matrix not found"
	SetDrawLayer UserBack
	TitleBox errorMessage,pos={17,28},size={344,16},title="Matrix not found. Remove cat from keyboard and try again."
	TitleBox errorMessage,fSize=14,frame=0
	Button buttonOk,pos={16,57},size={50,20},title="Ok",proc=okButton
	SetWindow nameError hook(xButton) = exitButtonHook
	PauseforUser nameError
else
	// Runs the decomposition procedure with the entered matrix name
	decomp(selectMatrix)
	
endif

end



//--------------------------------------------------------------------------------------------------------------
// Button control to handle the Compress button
//--------------------------------------------------------------------------------------------------------------

Function factCompressButton(ctrlName) :ButtonControl
String ctrlname

// Variable declaration
SVAR selectMatrix = root:selectMatrix, numbFact = root:numbFact
String dataFolder = "root:Factor_Analysis:" + selectMatrix +"_Factor_Analysis", selectedEigenValues = ""
Variable holder

// Error handling for factor compression. Checks is decomposition folder exists.
if(DataFolderExists(dataFolder) == 0)
	PauseUpdate; Silent 1		// building window...
	NewPanel /N=decompError/W=(150,77,556,167) as "Decomposition Error"
	TitleBox errorMessage,pos={17,28},size={344,16},title="Matrix doesn't exist or decomposition has not been performed."
	TitleBox errorMessage,fSize=14,frame=0
	Button buttonOk,pos={16,57},size={50,20},title="Ok",proc=okButton
	SetWindow decompError hook(xButton) = exitButtonHook
	PauseforUser decompError
	
else
	// Sets data folder to the decomposition data folder if it exists
	Setdatafolder $dataFolder
	
	// Variable and Wave declaration
	Variable n=0
	WAVE realEigenValues, eigenValueList = root:eigenValueList
	
	// Do-while loop to pick out the selected eigen values
	Do 
		holder = eigenValueList[n]
	
		If(holder >= 1)
			selectedEigenValues = selectedEigenValues + num2str(n) + ";"

		else
		
		endif

		n+=1
	
	While(n<numpnts(eigenValueList))
	
	// Error handling to make sure at least one eigen value was selected
	if(ItemsinList(selectedEigenValues) == 0)
		setdatafolder root:Factor_Analysis
		decomp(selectMatrix)

		PauseUpdate; Silent 1		// building window...
		NewPanel /N=eigenError/W=(150,77,416,167) as "Eigen Value Error"
		TitleBox errorMessage,pos={17,28},size={344,16},title="Please select at least one eigen value."
		TitleBox errorMessage,fSize=14,frame=0
		Button buttonOk,pos={16,57},size={50,20},title="Ok",proc=okButton
		SetWindow eigenError hook(xButton) = exitButtonHook
		PauseforUser eigenError
		
	else
		// Passes selected Eigen values and Matrix to the factor compression procedure
		numbFact = selectedEigenValues
		factCompress(selectMatrix,numbFact)
	
	endif

endif

end



//--------------------------------------------------------------------------------------------------------------
// Button Control Procedure to kill the various error windows when the user 
// presses the ok button.
//--------------------------------------------------------------------------------------------------------------

Function okButton(ctrlName) :ButtonControl
String ctrlName

DoWindow/C killError
Killwindow killError

end



//--------------------------------------------------------------------------------------------------------------
// Set variable control procedure to make sure Eigen Value list always matches
// the eigen values of the selected matrix if they have been calculated. If not
// it clears the eigen value list to prevent a user error.
//--------------------------------------------------------------------------------------------------------------

Function setVarCntrl(SV_Struct) :SetVariableControl
STRUCT WMSetVariableAction &SV_Struct

// Variable declaration
SVAR selectMatrix = root:selectMatrix
String dataFolder = "root:Factor_Analysis:" + selectMatrix +"_Factor_Analysis"

// Detects events that shouldn't update the eigen value list 
// and prevents the update from happening
if(SV_Struct.eventCode == -1 || SV_Struct.eventCode == 6)

else
	// Sets the data folder so the decomp procedure can be called
	Setdatafolder root:Factor_Analysis

	// Clears the list if the eigen values have not been calculated
	if(datafolderExists(dataFolder) == 0)
		Make/O/T/N=(0,5) root:eigenListBox
		Make/O/B/N=0 root:eigenValueList=0

	else
		// Populates the list with the correct eigen values if they have been calculated
		decomp(selectMatrix)

	endif

endif

end



//--------------------------------------------------------------------------------------------------------------
// Hook function to make the exit button in the upper right hand corner work
// like the cancel button
//--------------------------------------------------------------------------------------------------------------
Function exitButtonHook(s)
STRUCT WMWinHookStruct &s
Variable statusCode = 0

strswitch( s.eventName )

	case "killVote":
		strswitch( s.winName )
		
		case "faMainPanel":
			cancelButton("mainCancel")
			break
		
		case "dSelect":
			cancelButton("cancel")
			break
		
		endswitch
	
	if(strsearch(s.winName, "Error", 0) != -1)
		okButton("")
	endif	
		
	statusCode = 1	
	break
	
endswitch
return statusCode
end

//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//
//																	 //
//																	 //
//  						Main Procedures								 //
//																	 //
//																	 //
//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//


//--------------------------------------------------------------------------------------------------------------
// Function for creating fake data to test the procedures
//--------------------------------------------------------------------------------------------------------------

Function mockData(fwhm, number)
// Variables for various peak parameters
Variable fwhm, number
Variable sigma = fwhm/(2*sqrt(2*ln(2))),p1,p2, mu

NewDataFolder/O root:mockData
SetDataFolder root:mockData

// Makes data waves
Do
	String wName = "data" + num2str(p1+1)
	Make/O/N=1201 $wName
	Wave w3 = $wName
	mu = enoise(500)
	w3 = exp(-((x-(mu+dimsize(w3,0)/2))^2)/(2*(sigma^2)))/(sigma*sqrt(2*pi))
	//SetScale x -1,1, w3
	WaveStats/Q w3
	w3 = (w3-V_min)/(V_max-V_min)

	//sigma = (fwhm+fwhm*(enoise(1)))/(2*sqrt(2*log(2)))
	p1+=1
while(p1 < number)


// Creates a matrix of linear combinations of the data waves
variable i = 0
string setName


do
	setName = "set" + num2str(i+1)
	Make/O/N=1201 $setName = 0
	
	WAVE w1 = $setName
	p2 = 0
	Do
		String dataname = "data" + num2str(p2+1)
		Wave w2 = $dataname
		Variable x = enoise(.5)+.5
		w1 = w1 + x*w2
		
		p2+=1
	While(p2 < number)
	i += 1
	//WaveStats/Q w1
	//w1 = (w1-V_min)/(V_max-V_min)
	//w1 = enoise(.05)+w1
while (i<= number*10)

SetDataFolder root:

end



//--------------------------------------------------------------------------------------------------------------
// This procedure performs the Factor Decomposition of a selected Matrix
//--------------------------------------------------------------------------------------------------------------

Function deComp(Matrix)
// Variable declaration
String Matrix
NVAR normalizeDat = root:normalizeDat
Variable/G isTransposed = 0
// Creates the name of the data folder to hold the results
String dataFolderName = Matrix + "_Factor_Analysis"

// Gets the current data folder and stores the path to the martix
String savDF = Getdatafolder(1)
String pathToWave = savDF + Matrix

// Checks for the data folder associated with the input matrix and skips the calculation if it exists.
if(datafolderExists(dataFolderName))
	
	setdataFolder $dataFolderName
	
else

	// Creates a data folder to hold the results of the Factor Analysis
	NewDataFolder/O/S $dataFolderName
	
endif

	// Creates the Raw Data Martix that will be used in the calculations
	if(exists("rawData"))
	KillWaves cols,rows,rawData
	endif
	Wave dataMatrix = $pathToWave
	If(dimsize(dataMatrix,1)>dimsize(dataMatrix,0))
		MatrixTranspose dataMatrix
		isTransposed = 1
	endIf
	Wave rawData
	Duplicate dataMatrix, rawData 
	
	// Retrieves the numbuer of colums and rows in the dataMatrix
	MatrixOp/O cols=numCols(dataMatrix)
	MatrixOp/O rows=numRows(dataMatrix)

	// Calculates the Covariance Matrix
	MatrixOp/O covData = dataMatrix^t x dataMatrix
	//MatrixOp/O covData = subtractMean(covData,1)
	if(normalizeDat == 1)
		MatrixOp/O covData = NormalizeCols(covData)
	else
	endif

	// Calculates the Eigen Vectors of the Covariance Matrix
	MatrixEigenV/R covData

	// The Column factors are the Eigen Vectors of the Covariance Matrix
	Duplicate/O M_R_eigenVectors, colVect

	//Calculates the Row factors
	MatrixOp/O rowVect = dataMatrix x colVect


// Retrieves the Real Part of the eigen Values for display in the main control panel List Box
WAVE W_eigenValues
Make/O/N=(numpnts(W_eigenValues)) realEigenValues = real(W_eigenValues)
//Sort/R realEigenValues realEigenValues
Make/O/T/N=(numpnts(realeigenValues),6) root:eigenListBox
WAVE/T w12 = root:eigenListBox

MatrixOp/O ecols=numCols(realEigenValues)
MatrixOp/O erows=numRows(realEigenValues)
Variable numCols = ecols[0], numRows = erows[0]


//Creates Wave to hold the selection status of the main control panel List Box
Make/O/B/N=(numpnts(realeigenValues)) root:eigenValueList=0, xValues

//Creates Waves to hold the Error results
Make/O/N=(numpnts(realeigenValues)) REvalues, IEvalues, XEvalues, INDvalues, REVfvalues

//Calculates the various error functions for 1 through n factors
Variable n = 1, j = 0,  eigenValueHolder
String RError, IError, XError, IND, eigenValueString, REVf

Do
REvalues[n-1] = realError(Matrix, n)
IEvalues[n-1] = imbedError(Matrix, n)
XEvalues[n-1] = extractError(Matrix, n)
INDvalues[n-1] = factIND(Matrix, n)
REVfvalues[n-1] = redeigen(Matrix, n)
xValues[n-1] = n

eigenValueHolder = realeigenValues[n-1]
eigenValueString = num2str(eigenValueHolder)
RError = num2str(realError(Matrix, n))
IError = num2str(imbedError(Matrix, n))
XError = num2str(extractError(Matrix, n))
IND = num2str(factIND(Matrix,n))
REVf = num2str(redeigen(Matrix,n))
w12[n-1][0] = eigenValueString
w12[n-1][1] = RError
w12[n-1][2] = IError
w12[n-1][3] =  XError
w12[n-1][4] =  IND
w12[n-1][5] = REVf
n+=1
While((n-1) < numpnts(realeigenValues))

end



//--------------------------------------------------------------------------------------------------------------
// This procedure performs Factor compression of the selected matrix using 
// the selected Eigen Values
//--------------------------------------------------------------------------------------------------------------

Function factCompress(Matrix,numFact)
// Variable declaration 
String Matrix, numFact

// Stores the path to the Matrix to be compressed
String pathToMatrix = "root:Factor_Analysis:" + Matrix

// Declares the waves to be used in the Calculation. These were created during decomposition.
Wave dataMatrix =$pathToMatrix, W_eigenValues,colRed,rowVect, colVect, rawdata, datared

// Retrieves the number of columns and rows from the matrix and prints them
Variable numCols = dimsize(dataMatrix,1), numRows = dimsize(dataMatrix,0)

// Retrieves the real part of the eigen values
MatrixOp/O lamda = real(W_eigenValues)

// Reproduces the data using the entire factor space from the decomposition step
MatrixOp/O repData = rowVect x colVect^t

// Variable declaration of do-while loop that extracts the choosen column and row vectors
variable n = 0

// Creates the martices to hold the choosen column and row vectors
Make/O/N=(numCols,ItemsInList(numFact)) colRed
Make/O/N=(numRows,ItemsInList(numFact)) rowRed

// Variable decelartion for nested do-while loop to create column and row matrices
variable i=0

// do-while loop for creating the reduced column matrix
Do
	String setName = "colRed" + num2str(str2num(stringfromlist(i,numFact))+1)
	Wave w3= $setName
	colRed[][i] = colVect[p][str2num(stringfromlist(i,numFact))]
	
	i+=1
while (i < ItemsInList(numFact))

// Variable reinitilization
i = 0 

// do-while loop for creating the reduced row matrix
Do
	setName = "rowRed" + num2str(str2num(stringfromlist(i,numFact))+1)
	Wave w4 = $setName	
	rowRed[][i] = rowVect[p][str2num(stringfromlist(i,numFact))]

	i+=1
while ( i < ItemsInList(numFact))

// Reproduces the Data with the reduced Factor Space
MatrixOp/O dataRed = rowRed x colRed^t

Variable p1 = 0
String dimLabel
Do
If(stringMatch(GetDimLabel(dataMatrix,1,p1),""))
	dimLabel = "datared" + num2str(p1+1)
else
	dimLabel = GetDimLabel(dataMatrix,1,p1)[0,26] + "_Red"	
endif
SetDimLabel 1,p1,$dimLabel,dataRed
p1+=1
while(p1 < dimsize(dataRed,1))

MatrixOp/O dataDiff = dataMatrix - dataRed

// Retrieves the number of rows for window positioning
Variable winWidth = dimSize(dataMatrix,0) * 1.25

// Displays the raw data
NewImage/K=1 rawdata
ModifyImage rawdata ctab={*,*,Rainbow,0}
//ModifyGraph width = {perUnit, 1.25, top}, height = {Aspect, 1}
MoveWindow 0,0,400,400

// Displays the reproduced data using the reduced Factor Space
NewImage/K=1 dataRed
ModifyImage dataRed ctab={*,*,Rainbow,0}
//ModifyGraph width = {perUnit, 1.25, top}, height = {Aspect, 1}
MoveWindow 450,0,850,400

splitMat("rawData")
splitMat("dataRed")

String dataFolderName = "factors"
Variable p2 = 0
Do
dataFolderName = dataFolderName + "_" + num2str(str2num(stringfromlist(p2,numFact))+1)
p2 += 1
while(p2 < ItemsInList(numFact))

NewDataFolder/O/S $dataFolderName
Duplicate/O ::rawData,:rawData
Duplicate/O ::dataRed,:dataRed
Duplicate/O ::dataDiff,:dataDiff

splitMat("rawData")
splitMat("dataRed")
SetDataFolder ::
end

//--------------------------------------------------------------------------------------------------------------
// Procedure for target testing the data.
//--------------------------------------------------------------------------------------------------------------

Function targetTest(numFact, wave0)
// Variable Declaration
variable numFact
Wave wave0

// Sets the matrix to be used in the calcuations
Wave rawData

// Gets the number of columns and rows of the Data matrix
MatrixOp/O cols=numCols(rawData)
//MatrixOp/O rows=numRows(rawData)
Variable numCols = cols[0], numRows = numpnts(wave0)
KillWaves cols,rows

// Makes the wave to hold the test data matrix
Make/O/N=(numRows) testVect

// Declares the waves to be used for target testing
Wave data1,data2, rowRed, colRed, rowRed1, colRed1

// Variable Declaration of do-while loop to create test matrix
Variable j=0

// do-while loop to create test matrix
do
	//if(j == numfact)
		testVect[j] = wave0[j]
	//else
	//	testData[0][j] = 0
	//endif
	//testData[j][1] = wave2[j]
	j+=1
while((j+1)<=numRows)


//Calculates the Transformation Matrix
MatrixOP/O trans = (Inv(rowRed^t x rowRed) x rowRed^t) x testVect
//MatrixOP/O trans = testData x colRed

//Calculates the Real Factors from the transformation matrix
//MatrixOp/O xData = rowRed x Inv(trans)
//MatrixOp/O yData = Inv(trans) x colRed^t

//Reproduces the Data with the Real Factors
//MatrixOp/O targetData = xData x yData
MatrixOp/O targetVect = rowRed x trans

end

//--------------------------------------------------------------------------------------------------------------
// Target testing interface
//--------------------------------------------------------------------------------------------------------------
Function testInt(numFact, targetList)
Variable numFact
String targetList

Wave rowRed, colRed

Make/T/O/N=(ItemsInList(targetList),6) ttestingRes

Variable n = ItemsInList(targetList), j = 0, k = 0
String stringHolder
Do
stringHolder = stringFromList(j, targetList)
//String nWave = rescale(stringHolder, rawdata1)
String nWave = stringHolder
targetTest(numFact, $nWave)

ttestingRes[j][0] = nWave
ttestingRes[j][1] = num2str(apparentError(numFact))
ttestingRes[j][2] = num2str(rootMS(numFact))
ttestingRes[j][3] = num2str(realErrorTest(numFact))
ttestingRes[j][4] = num2str(SPOIL(numFact))
ttestingRes[j][5] = num2str(testSL(numFact))

Make/O/N=(numpnts($nWave), numfact) testData

k = 0
Wave testVect
Do
testData[k][j] = testVect[k]

k+=1
while(k < numpnts($nWave))

j+=1
while(j < n)

MatrixOP/O transF = (Inv(rowRed^t x rowRed) x rowRed^t) x testData

MatrixOp/O xData = rowRed x transF
MatrixOp/O yData = Inv(transF) x colRed^t
//xData = abs(xData)

MatrixOp/O targetData = xData x yData

Wave targetData
String nameHolder2
MatrixOp/O cols=numCols(targetData)
MatrixOp/O rows=numRows(targetData)
Variable numCols = cols[0], numRows = rows[0]

Variable l = 0
Do
	Variable q = 0
	nameHolder2 = "targetdata" + num2str(l+1)
	Make/O/N=(numrows) $nameHolder2
	Wave w1 = $nameHolder2
	Do
		w1[q] = targetData[q][l]
		q+=1
	while(q<=numRows)
l+=1
while(l<numcols)

end

//--------------------------------------------------------------------------------------------------------------
// Not sure what this is supposed to do. Look into that.
//--------------------------------------------------------------------------------------------------------------

Function keyset()
// Normalizes the row factor space for some reason.
Wave rowVect
Matrixop rowVect_N = NormalizeCols(rowVect)

end



//--------------------------------------------------------------------------------------------------------------
// Procedure for creating the data matrix for factor analysis.
//--------------------------------------------------------------------------------------------------------------

Function dataMatrix(wList, matrixName)
// Variable Declaration
String wList, matrixName

// Creates the Folder to hold all the factor Anaylsis stuff
NewDataFolder/S/O Factor_Analysis

// Variable Initilization for do-while loop to create the data matrix
Variable j = itemsInList(wList), i = 0, k = 0

// Makes the empty data Matrix Wave
WaveStats/Q $stringFromlist(0, wList)
Make/O/N=(V_npnts, j) $matrixName
Wave w1 = $matrixName

//Do while loop to put the data in the empty data Matrix
Do
	String dimLabel = stringFromList(k,wList)[strsearch(stringFromList(k,wList),":",strsearch(stringFromList(k,wList),":",inf,1)-1,1)+1,inf]
	dimLabel = replacestring("_", dimLabel, "")
	dimLabel = replacestring(":", dimLabel, "_")
	dimLabel = replacestring("'", dimLabel, "")
	dimLabel = replacestring("txt",dimLabel,"")
	dimLabel = CleanupName(dimLabel, 0)
	SetDimLabel 1,k,$dimLabel,w1
	Wave w2 = $stringFromList(k, wList)
	Setscale x,pnt2x(w2,0),pnt2x(w2,numpnts(w2)),w1
	w1[][k] = w2[p]
	k+=1
while(k<j)

//Old code I probably don't need but I'll leave it here for now
//--------------------------------------------------------------------------------------------------------------
//String energy = stringFromlist(1, wList)
//Variable repString = strsearch(energy,":", inf, 1)
//energy = replacestring(energy[repstring+1, inf],energy,"cor_E")

//Wave w3 = $energy
//WaveStats/Q w3
//Setscale x V_min,V_max, w1
//----------------------------------------------------------------------------------------------------------

// Creates matrix for f(z) color scale when waterfall plotting
String dataMatrix_y = matrixName + "_y"
Duplicate/O w1, $dataMatrix_y

Wave w4 = $dataMatrix_y
w4 = y

// Plots the data matrix
NewWaterFall/K=1 w1
ModifyGraph zColor($matrixName)={$dataMatrix_y,*,*,Rainbow,0}

end


//--------------------------------------------------------------------------------------------------------------
// Procedure for testing for uniqueness and constant factors
//--------------------------------------------------------------------------------------------------------------

Function ucTest()

// Sets the matrix to be used in the calcuations
Wave rawData

// Gets the number of columns and rows of the Data matrix
MatrixOp/O cols=numCols(rawData)
MatrixOp/O rows=numRows(rawData)
Variable numCols = cols[0], numRows = rows[0]
KillWaves cols,rows

// Makes the wave to hold the test data matrix
Make/O/N=(numRows, 1) utestCol
Make/O/N=(1, numCols) utestRow
Make/O/N=(numRows, 1) ctestCol
Make/O/N=(1, numCols) ctestRow

// Declares the waves to be used for target testing
Wave data1,data2, rowRed, colRed, rowRed1, colRed1, wave0, wave2

// Variable Declaration of do-while loop to create test matrix
Variable j=0, k = 0

// do-while loop to create test matrix
do
	ctestCol[j][0] = 1
	j+=1
	
while((j+1)<=numRows)

do
	ctestRow[0][k] = 1
	k+=1
	
while((k+1)<=numCols)

//Calculates the Transformation Matrix
MatrixOP/O ctransCol = (Inv(rowRed^t x rowRed) x rowRed^t) x ctestCol
MatrixOP/O ctransRow = ctestRow x colRed

//Reproduces the Data with the Real Factors
MatrixOp/O ctargetCol = rowRed x ctransCol
MatrixOp/O ctargetRow = ctransRow x colRed^t

Variable l=0, m=0, urow=0, ucol=0
String saveDF = GetDataFolder(1)

NewDataFolder/O/S utestrow
do
	String rownameHolder = "utestRow_" + num2str(urow+1)
	String rowdataName = "utargetRow_" + num2str(urow+1)

	Make/O/N=(numRows,1) $rownameHolder
	Make/O/N=(numRows,1) $rowdataName
	Wave w5 = $rownameHolder, w6 = $rowdataName
	l = 0
	// do-while loop to create test matrix
	do
		if(l == urow)
			w5[l][0] = 1
		else
			w5[l][0] = 0
		endif
		
		l+=1
	
	while((l+1)<=numRows)
	
	//Calculates the Transformation Matrix
	MatrixOP/O ctransCol = (Inv(rowRed^t x rowRed) x rowRed^t) x w5

	//Reproduces the Data with the Real Factors
	MatrixOp/O w6 = rowRed x ctransCol


	urow+=1

	if((urow+1) > numRows)
		SetDataFolder saveDF
	endif
	
while((urow+1)<=numRows)

saveDF = GetDataFolder(1)

NewDataFolder/O/S utestCol
do
	String colnameHolder = "utestCol_" + num2str(ucol+1)
	String coldataName = "utargetCol_" + num2str(ucol+1)

	Make/O/N=(1, numCols) $colnameHolder
	Make/O/N=(1, numCols) $coldataName
	Wave w7 = $colnameHolder, w8 = $colDataName
	m = 0
	
// do-while loop to create test matrix
	do
		if(m == ucol)
			w7[0][m] = 1
		else
			w7[0][m] = 0
		endif
		
		m+=1
	
	while((m+1)<=numCols)
	
	//Calculates the Transformation Matrix
	MatrixOP/O ctransRow = w7 x colRed

	//Reproduces the Data with the Real Factors
	MatrixOp/O w8 = ctransRow x colRed^t


	ucol+=1

	if((ucol+1) > numCols)
		SetDataFolder saveDF
	endif
	
while((ucol+1)<=numCols)

end

//Procedure for creating a random data matrix
Function randMatrix(matrixName, scale)
String matrixName
Variable scale
WAVE w1 = $matrixName

// Gets the number of columns and rows of the Data matrix
MatrixOp/O cols=numCols(w1)
MatrixOp/O rows=numRows(w1)
Variable numCols = cols[0], numRows = rows[0]
KillWaves cols,rows
Make/O/N=(numRows,numCols) randomMatrix
Variable j=0,k=0

Do
	Do
		randomMatrix[j][k] = gnoise(1)/scale
		j+=1
	while(j < numRows)
	j = 0
	k+=1
while(k < numCols)

end

//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//
//																	 //
//																	 //
//  					Error Testing Procedures							 //
//																	 //
//																	 //
//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//

//--------------------------------------------------------------------------------------------------------------
// Procedure for determining Real Error (RE or RSD)
//--------------------------------------------------------------------------------------------------------------
Function realError(matrixName, n)
String matrixName
Variable n

Wave w10 = $matrixName, realEigenValues, rawData
MatrixOp/O cols=numCols(rawData)
MatrixOp/O rows=numRows(rawData)
Variable numCols = cols[0], numRows = rows[0], s, l

if(numCols > numRows)
	s = numRows
	l = numCols
elseif(numRows >= numCols)
	s = numCols
	l = numRows
endif

Variable sumEV, RError
sumEV = sum(realEigenValues,pnt2x(realEigenValues, n),pnt2x(realEigenValues, s))

RError = (sumEV/(l*(s-n)))^(1/2)

Return RError
end

//--------------------------------------------------------------------------------------------------------------
// Procedure for determining Imbeded Error (IE)
//--------------------------------------------------------------------------------------------------------------
Function ImbedError(matrixName, n)
String matrixName
Variable n

Wave w10 = $matrixName, rawData
MatrixOp/O cols=numCols(rawData)
MatrixOp/O rows=numRows(rawData)
Variable numCols = cols[0], numRows = rows[0], s, l

if(numCols > numRows)
	s = numRows
	l = numCols
elseif(numRows >= numCols)
	s = numCols
	l = numRows
endif

Variable IError

IError = realError(matrixName, n)*sqrt(n/s)

Return IError
end

//--------------------------------------------------------------------------------------------------------------
// Procedure for determining Extracted Error (XE)
//--------------------------------------------------------------------------------------------------------------
Function extractError(matrixName, n)
String matrixName
Variable n

Wave w10 = $matrixName, rawData
MatrixOp/O cols=numCols(rawData)
MatrixOp/O rows=numRows(rawData)
Variable numCols = cols[0], numRows = rows[0], s, l

if(numCols > numRows)
	s = numRows
	l = numCols
elseif(numRows >= numCols)
	s = numCols
	l = numRows
endif

Variable XError

XError = realError(matrixName, n)*sqrt((s-n)/s)

Return XError
end

//--------------------------------------------------------------------------------------------------------------
// Procedure for determining Factor Indicator Function (IND)
//--------------------------------------------------------------------------------------------------------------
Function factInd(matrixName, n)
String matrixName
Variable n

Wave w10 = $matrixName, rawData
MatrixOp/O cols=numCols(rawData)
MatrixOp/O rows=numRows(rawData)
Variable numCols = cols[0], numRows = rows[0], s, l

if(numCols > numRows)
	s = numRows
	l = numCols
elseif(numRows >= numCols)
	s = numCols
	l = numRows
endif

Variable IND

IND = realError(matrixName, n)/(s-n)^2

Return IND
end

//--------------------------------------------------------------------------------------------------------------
//Procedure for determing Reduced Eigenvalues and significance level
//--------------------------------------------------------------------------------------------------------------
Function redeigen(matrixName,n)
String matrixName
Variable n

Wave w10 = $matrixName, rawData, realEigenValues
Make/O/N=(numpnts(realEigenValues)) revEigenValues
MatrixOp/O cols=numCols(rawData)
MatrixOp/O rows=numRows(rawData)
Variable numCols = cols[0], numRows = rows[0], s, l

if(numCols > numRows)
	s = numRows
	l = numCols
elseif(numRows >= numCols)
	s = numCols
	l = numRows
endif

Variable j=0, revN, Fvalue

Do
revEigenValues[j] = realEigenValues[j]/((numRows-(j+1)+1)*(numcols-(j+1)+1))
j+=1
while(j<numpnts(realEigenValues))

revN = revEigenValues[n-1]

Fvalue = revN*(s-n)/sum(revEigenValues, n, s)
Print "F: ",Fvalue,"s-n: ",s-n,"FCDF: ",StatsFCDF(Fvalue, 1,s - n)
Return (1-StatsFCDF(Fvalue,1,s - n))*100
Return (1 - ((s-n)*Fvalue/((s-n)*Fvalue +1)))*100
end

//--------------------------------------------------------------------------------------------------------------
// Procedure for determining Apparent Error in the Test Vector (AET)
//--------------------------------------------------------------------------------------------------------------
Function apparentError(numFacts)
Variable numFacts
Wave testVect, targetVect
Make/O/N=(numpnts(targetVect)) difference
difference = targetVect - testVect

Return ((sum(difference))^2/(numpnts(targetVect) - numfacts))^(.5)

end

//--------------------------------------------------------------------------------------------------------------
// Procedure for determining Root Mean Square of the Predicted Vector (REP)
//--------------------------------------------------------------------------------------------------------------
Function rootMS(numFacts)
Variable numFacts
Wave REvalues, trans

Return REvalues[numFacts-1]*norm(trans)
end

//--------------------------------------------------------------------------------------------------------------
// Procedure for determining Real Error in the Traget Vector (RET)
//--------------------------------------------------------------------------------------------------------------
Function realErrorTest(numFacts)
Variable numFacts

Return ((apparentError(numfacts))^2-(rootMS(numfacts))^2)^(.5)
end

//--------------------------------------------------------------------------------------------------------------
// Procedure for determining the SPOIL function
//--------------------------------------------------------------------------------------------------------------
Function SPOIL(numFacts)
Variable numFacts

Return realErrorTest(numFacts)/(rootMS(numFacts))
end 

//--------------------------------------------------------------------------------------------------------------
// Procedure for determining the percent significance level of the test vector
//--------------------------------------------------------------------------------------------------------------
Function testSL(numfacts)
Variable numFacts
Wave rawData
MatrixOp/O cols=numCols(rawData)
MatrixOp/O rows=numRows(rawData)
Variable numCols = cols[0], numRows = rows[0], s, l

if(numCols > numRows)
	s = numRows
	l = numCols
elseif(numRows >= numCols)
	s = numCols
	l = numRows
endif

Variable j = numFacts+1, sum01 = 0, Fvalue
Do
sum01 = sum01 + (numrows - j + 1)*(numcols-j+1)

j+=1
while(j<=s)

Fvalue = ((apparentError(numFacts)/rootMS(numfacts))^(2))*(sum01/((numrows-numfacts+1)*(numcols-numfacts+1)))

Return (1-StatsFCDF(Fvalue,l - numFacts, s - numFacts))*100
Return (1-((s-numfacts)*Fvalue/((s-numfacts)*Fvalue +(numcols-numfacts))))*100
end


//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//
//																	 //
//																	 //
//  						    Utility Procedures							 //
//																	 //
//																	 //
//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//


//--------------------------------------------------------------------------------------------------------------
// Procedure for deconstructing the data matrices
//--------------------------------------------------------------------------------------------------------------
Function splitMat(matName)
String matName
String nameHolder1, nameHolder2
Wave w1 = $matName

Variable k = 0

Do
	nameHolder1 = matName + "_" + num2str(k+1)
	if(stringmatch(GetDimLabel(w1,1,k), "") == 0)
		nameHolder2 = GetDimLabel(w1,1,k)
		Make/O/N=(dimsize(w1,0)) $nameHolder2
		Wave w3 = $nameHolder2
		Setscale/P x,dimoffset(w1,0),dimdelta(w1,0), w3
		w3[] = w1[p][k]
	endif
	Make/O/N=(dimsize(w1,0)) $nameHolder1
	Wave w2 = $nameHolder1
	Setscale/P x,dimoffset(w1,0),dimdelta(w1,0), w2
	w2[] = w1[p][k]
	k+=1
while(k<dimsize(w1,1))

end

Function splitMatCols(matName)
String matName
Wave w1 = $matName
Make/O/N=(dimsize(w1,0),dimsize(w1,1)) $(matName + "_trans")
Wave temp = $(matName + "_trans")

MatrixOp/O temp = w1^t

String tempName = NameOfWave(temp)

splitmat(tempName)

end

//--------------------------------------------------------------------------------------------------------------
// Procedure for rescalling test vectors
//--------------------------------------------------------------------------------------------------------------
Function/S rescale(testVectS, rawData1)
String testVectS
Wave rawData1

Wave testVect = $testVectS
String nWave = testVectS[strsearch(testVectS, ":", inf, 1)+1,inf]
//nWave = CleanupName(nWave, 0)

KillWaves/Z $nWave
Duplicate testVect $nWave
Wave testVect = $nWave
SetScale x, 0, (pnt2x(testVect, numpnts(testVect)-1) - pnt2x(testVect, 0)), testVect
Resample/SAME=rawData1 testVect



Wavestats/Q testVect
Variable maxRowLoc = V_maxRowLoc, maxRowLocT = trunc(pnt2x(testVect, maxRowLoc)), maxTest = V_max
WaveStats/Q rawData1
Variable maxRaw = V_max

InsertPoints 0, (maxRowLocT-maxRowLoc), testVect

Setscale x, 0, numpnts(rawData1), testVect
Redimension/N=(numpnts(rawData1)) testVect

Return nWave
end

//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//
//																	 //
//																	 //
//  				        Self Modeling Curve Analysis						 //
//																	 //
//																	 //
//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//

Function smcAnalysis()
Wave dataRed
Make/O/N=(dimsize(dataRed,0),dimsize(dataRed,1)) secondMoment
Variable n = dimsize(dataRed,0)

//if(dimsize(dataRed,0) > 10000)
//Duplicate/O dataRed, dataRed_samp
//Resample/DOWN=10 dataRed_samp
//MatrixOp/O secondMoment = (dataRed_samp x dataRed_samp^t)/n

//if(dimsize(rowRed,0) > 10000)
//Resample/DOWN=10 rowRed
//endif
//else
MatrixOp/O secondMoment = (dataRed x dataRed^t)/n
MatrixEigenV/L secondMoment
//MatrixEigenV/SYM/EVEC/RNG={2,(dimsize(secondMoment,0)-1),(dimsize(secondMoment,0))} secondMoment
//endif

Wave M_eigenVectors,M_L_eigenVectors

Variable i = 0
Make/O/N=(dimsize(dataRed,1)) eta1,eta2

Make/O/N=(dimsize(secondMoment,0)) eigen1,eigen2
Variable k = 0
Do 
//eigen1[k] = -M_eigenVectors[k][dimsize(M_eigenVectors,1)-1]
//eigen2[k] = -M_eigenVectors[k][dimsize(M_eigenVectors,1)-2]
if(M_L_eigenVectors[0][0] <= 0)
eigen1[k] = -M_L_eigenVectors[k][0]
eigen2[k] = -M_L_eigenVectors[k][1]
else
eigen1[k] = M_L_eigenVectors[k][0]
eigen2[k] = M_L_eigenVectors[k][1]
endif
k+=1
while(k < dimsize(secondMoment, 0))

Do
String nameHolder = "datared" + num2str(i+1)
//if(dimsize(dataRed,0) > 10000)
//String nameHolder1 = nameHolder + "_samp"
//Duplicate/O $nameHolder, $nameHolder1
//Resample/DOWN=10 $nameHolder1
//nameHolder = nameHolder1 
//endif
Wave w1 = $nameHolder
eta1[i] = MatrixDot(eigen1, w1)
eta2[i] = MatrixDot(eigen2, w1)

i+=1
while(i < dimsize(dataRed,1))

Display/B=center eta2 vs eta1
ModifyGraph mode=3,marker=19,msize=2
ModifyGraph standoff=0,freePos(center)={0,left}
ModifyGraph tick=1
Wavestats/Q eta2

If(V_min < 0)
SetAxis left (V_min + V_min/1.5),(V_max + V_max/1.5)
else
SetAxis left (V_min - V_min/1.5),(V_max + V_max/1.5)
endif

Wavestats/Q eta1
SetAxis center 0,(V_max + V_max/2)

Make/O/N=0 zetaWave,tauWave

k = 0
Variable neg = 0, pos = 0
Do
If(eigen2[k] >= 0)
InsertPoints pos+1, 1, zetaWave
zetaWave[pos] = abs(eigen1[k]/eigen2[k])
pos+=1
else
InsertPoints neg+1, 1, tauWave
tauWave[neg] = abs(eigen1[k]/eigen2[k])
neg+=1
endif
k+=1
while(k < numpnts(eigen2))

Variable zeta, tau

Wavestats/Q zetaWave
zeta = -V_min

Wavestats/Q tauWave
tau = V_min

Print zeta,tau

Variable c1,c2

c1 = sum(eigen1)
c2 = sum(eigen2)

Variable etaV1x = 1/(c1 + c2*zeta)
Variable etaV1y = zeta*etaV1x

Make/O/N=(numpnts(eta1)) etaRatio = eta2/eta1
Wavestats/Q etaRatio

Variable etaMax = V_max
Variable etaMin = V_min

Variable etaV2x = 1/(c1+c2*etaMax)
Variable etaV2y = etaMax*etaV2x

Variable etaV3x = 1/(c1+c2*etaMin)
Variable etaV3y = etaMin*etaV3x

Variable etaV4x = 1/(c1+c2*tau)
Variable etaV4y = tau*etaV4x

Make/O/N=(dimsize(secondMoment,0)) factor1Max, factor1Min, factor2Max, factor2Min

factor2Max = etaV1x*eigen1 + etaV1y*eigen2
factor1Max = etaV2x*eigen1 + etaV2y*eigen2
factor2Min = etaV3x*eigen1 + etaV3y*eigen2
factor1Min = etaV4x*eigen1 + etaV4y*eigen2

Wavestats/Q factor1Max
factor1Max = (factor1Max-V_min)/(V_max-V_min)
Wavestats/Q factor1Min
factor1Min = (factor1Min-V_min)/(V_max-V_min)
Wavestats/Q factor2Max
factor2Max = (factor2Max-V_min)/(V_max-V_min)
Wavestats/Q factor2Min
factor2Min = (factor2Min-V_min)/(V_max-V_min)

Display factor1Max, factor1Min
ModifyGraph mode(factor1Max)=7
ModifyGraph toMode(factor1Max)=1
ModifyGraph hbFill(factor1Max)=2
Display factor2Max, factor2Min
ModifyGraph mode(factor2Max)=7
ModifyGraph toMode(factor2Max)=1
ModifyGraph hbFill(factor2Max)=3
ModifyGraph usePlusRGB(factor2Max)=1,plusRGB(factor2Max)=(0,0,52224)
ModifyGraph useNegRGB(factor2Max)=1,negRGB(factor2Max)=(0,0,52224)
ModifyGraph rgb(factor2max)=(0,0,52224)
ModifyGraph rgb(factor2min)=(0,0,52224)

testInt(2,"factor1Max;factor2Max")
Wave ydata
Wave/T ttestingRes
Variable p1 =0

Make/O/N=(dimsize(secondMoment,0)) factor1test, factor2test
Make/O/N=10 eta1xtest, eta1ytest, eta2xtest, eta2ytest
Make/O/N=(dimsize(ydata,0),dimsize(ydata,1), dimsize(dataRed,1), 10) coefValues
Make/O/N=(dimsize(ttestingRes,0), dimsize(dataRed,1),10) resultValues

Variable eta1Step = (etaV2x - etaV4x)/10
Variable eta2Step = (etaV1x - etaV3x)/10

Do
eta1xtest[p1] = etaV2x + p1*eta1Step
eta1ytest[p1] = (1 - c1*eta1xtest[p1])/c2

eta2xtest[p1] = etaV1x + p1*eta2Step
eta2ytest[p1] = (1-c1*eta2xtest[p1])/c2

p1+=1
while(p1 < 10)

Variable p2 = 0, p3 = 0

Do
	factor1test = eta1xtest[p2]*eigen1 + eta1ytest[p2]*eigen2
	WaveStats/Q factor1test
	factor1test = (factor1test-V_min)/(V_max-V_min)
	p3 = 0
	Do
		factor2test = eta2xtest[p3]*eigen1 + eta2ytest[p3]*eigen2
		WaveStats/Q factor2test
		factor2test = (factor2test-V_min)/(V_max-V_min)
		testInt(2,"factor1test;factor2test")

		Variable p4 = 0
		Do
			Variable p5 = 0
			Do
				coefValues[p4][p5][p2][p3] = ydata[p4][p5]
				resultValues[p4][p2][p3] = str2num(ttestingRes[p4][5])
				p5+=1
			while(p5 < dimsize(ydata,1))
			p4+=1
		while(p4 < dimsize(ydata,0))

		p3+=1 
	while(p3 < 10)
p2+=1
while(p2 < dimsize(dataRed,1))

Make/O/N=(dimsize(dataRed,1)) scale1max, scale1min, scale2max, scale2min

Variable p6=0,p7=0,p8=0
Variable maxValue, minValue, difference, difference2
Do
	String nameHolder3, nameHolder4
	MatrixOp/O beamMat1 = beam(coefValues,0,p6)
	MatrixOp/O beamMat2 = beam(coefValues,1,p6)
	p7 = 0
	Do
		nameHolder3 = "scale1_" + num2str(p6) +"_" + num2str(p7)
		nameHolder4 = "scale1_" + num2str(p6)
		MatrixOp/O $nameHolder3 = row(beamMat1,p6)^t
		Wave wave1 = $nameHolder3, wave2 = $nameHolder4
		if(p7==0)
			//difference = abs(WaveMax($nameHolder3) - WaveMin($nameHolder3))
			difference = abs(99-resultValues[0][p6][p7])
			Make/O/N=(numpnts($nameHolder3)) $nameHolder4 = wave1
			KillWaves/Z $nameHolder3
		else
			//difference2 = abs(WaveMax($nameHolder3) - WaveMin($nameHolder3)) 
			difference2 = abs(99-resultValues[0][p6][p7]) 
			if(difference > difference2)
				Wave wave3 = $nameHolder3
				wave2 = wave3
				difference = difference2
				KillWaves/Z $nameHolder3
				nameHolder3 = "scale1_" + num2str(p6) +"_" + num2str(p7-1)
				KillWaves/Z  $nameHolder3
			elseif(difference < difference2)
				nameHolder3 = "scale1_" + num2str(p6) +"_" + num2str(p7)
				KillWaves/Z  $nameHolder3
			endif
		endif
		
		nameHolder3 = "scale2_" + num2str(p6) +"_" + num2str(p7)
		nameHolder4 = "scale2_" + num2str(p6)
		MatrixOp/O $nameHolder3 = row(beamMat2,p6)^t
		Wave wave1 = $nameHolder3, wave2 = $nameHolder4
		if(p7==0)
			//difference = abs(WaveMax($nameHolder3) - WaveMin($nameHolder3))
			difference = abs(99-resultValues[1][p6][p7])
			Make/O/N=(numpnts($nameHolder3)) $nameHolder4 = wave1
			KillWaves/Z $nameHolder3
		else
	 		//difference2 = abs(WaveMax($nameHolder3) - WaveMin($nameHolder3)) 
	 		difference2 = abs(99-resultValues[1][p6][p7]) 
			if(difference > difference2)
				Wave wave3 = $nameHolder3
				wave2 = wave3
				difference = difference2
				KillWaves/Z $nameHolder3
				nameHolder3 = "scale2_" + num2str(p6) +"_" + num2str(p7-1)
				KillWaves/Z  $nameHolder3
			elseif(difference < difference2)
				nameHolder3 = "scale2_" + num2str(p6) +"_" + num2str(p7)
				KillWaves/Z  $nameHolder3
			endif
		endif
		
		p7+=1
	while(p7 < 10)
	p6+=1
while(p6 < dimsize(dataRed,1))
Do
String nameHolder5 = "scale1_" + num2str(p8)
String nameHolder6 = "scale2_" + num2str(p8)
scale1max[p8] = WaveMax($nameHolder5)
scale1min[p8] = WaveMin($nameHolder5)
scale2max[p8] = WaveMax($nameHolder6)
scale2min[p8] = WaveMin($nameHolder6)
p8+=1
while(p8 < dimsize(datared,1))

Make/O/N=(numpnts(eigen1)) fact1,fact2

fact1 = eta1xtest[8]*eigen1 + eta1ytest[8]*eigen2
fact2 = eta1xtest[8]*eigen1 + eta1ytest[8]*eigen2

end

//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//
//																	 //
//																	 //
//		 Self Modeling Curve Resolution/2D Correlation Spectroscopy		 //
//																	 //
//																	 //
//--------------------------------------------------------------------------------------------------------------//
//--------------------------------------------------------------------------------------------------------------//

Function/S cs2d(input, high, low)
String input
Variable high,low
Variable p1,p2,p3
Wave inputMat = $input

Make/O/N=(dimsize(inputMat,0),dimsize(inputMat,1)) dynamicSpec, referenceSpec
Make/O/N=(dimsize(inputMat,1),dimsize(inputMat,1)) hntMat

Do

	Make/O/N=(dimsize(inputMat,1)) workVect
	workVect[p] = inputMat[p1][p]
	referenceSpec[p1][] = sum(workVect)/(high - low)

	p1+=1
while(p1 < dimsize(inputMat,0))

dynamicSpec = inputMat - referenceSpec

Variable m = dimsize(inputMat,1)

MatrixOp/O synCorr = (dynamicSpec x dynamicSpec^t)
synCorr = 1/(m -1) * synCorr

 Do
 	p3 = 0
 	Do
 		if(p2 == p3)
 			hntMat[p2][p3] = 0
 		else
 			hntMat[p2][p3] = 1/(pi*(p3 - p2))
 		endif
 		p3+=1
 	while(p3 < dimsize(inputMat,1))
 
 	p2+=1
 while(p2 < dimsize(inputMat,1))
 
MatrixOp/O asynCorr =(dynamicSpec x hntMat x dynamicSpec^t)

asynCorr = 1/(m - 1) * asynCorr

Variable peak1x, peak1y, peak2x, peak2y
FindValue/Z/V=(waveMax(asynCorr)) asynCorr

peak1x = floor(V_value/dimsize(asynCorr,0))
peak1y = V_value-peak1x*dimsize(asynCorr,0)

FindValue/Z/V=(waveMin(asynCorr)) asynCorr

peak2x =  floor(V_value/dimsize(asynCorr,0))
peak2y = V_value-peak2x*dimsize(asynCorr,0)

//String locations = num2str(peak1x) + ";" + num2str(peak1y)

String locations = "118,145,195,308"
//"91;142"
//;155;173;206;169,108,89"

Return locations

end

Function alsSMCR(input, locations, high, low, iterations)
String input, locations
Variable high,low,iterations
Wave inputMat = $input

//String locations = "118;145;195;308"

KillWaves/Z transMat,spectMat
Make/O/N=(dimsize(inputMat,1),itemsinlist(locations)) transMat

Variable p1,p2,p3,p4,p5,p6

Do
	transMat[][p1] = inputMat [str2num(stringfromlist(p1,locations))][p]
	p1+=1
while(p1 < itemsinlist(locations))

Do
	p2 = 0
	p3 = 0
	
	MatrixOp/O spectMat = inv(transMat^t x transMat) x transMat^t x inputMat^t
	MatrixOp/O spectMat = spectMat^t

	Do
		p3 = 0
		Do
			if(spectMat[p2][p3] < 0)
				spectMat[p2][p3] = 0
			endif
			
			p3+=1
		while(p3 < dimsize(spectMat,1))

		p2+=1
	while(p2 < dimsize(spectMat,0))

	MatrixOp/O transMat = inv(spectMat^t x spectMat) x spectMat^t x inputMat
	MatrixOp/O transMat = transMat^t
	MatrixOp/O repData = spectMat x transMat^t

	p4+=1
while(p4 < iterations)



Do
	Make/O/N=(dimsize(transMat,1)) workVect
	workVect[] = transMat[p5][p]
	p6 = 0
	Do
		transMat[p5][p6] = transMat[p5][p6]/sum(workVect)
	
		p6+=1
	while(p6 < dimsize(transMat,1))
	
	p5+=1
while(p5 < dimsize(transMat,0))

splitMat("spectMat")
splitMat("transMat")
end